R Notebook

# Lista de pacotes necessários
packages <- c(
  "dplyr", "ggplot2", "corrplot", "car", "lubridate",
  "tidyr", "stringr", "stringi", "rstatix", "sf",
  "purrr", "GGally", "FactoMineR", "factoextra",
  "patchwork", "mice", "janitor", "scales", "trend"
)

# Instalar apenas os que não estão instalados
installed <- packages %in% rownames(installed.packages())
if (any(!installed)) {
  install.packages(packages[!installed])
}

library(dplyr)
library(ggplot2)
library(corrplot)
library(car)
library(lubridate)
library(tidyr)
library(stringr)
library(stringi)
library(rstatix)
library(sf)
library(purrr)
library(GGally)
library(FactoMineR)
library(factoextra)
library(patchwork)
library(mice)
library(janitor)
library(scales)
library(trend)

#ler data
hypertension <- read.csv2("https://transparencia.sns.gov.pt/api/explore/v2.1/catalog/datasets/hipertensao/exports/csv?lang=pt&timezone=Europe%2FLisbon")
avc <- read.csv2("https://transparencia.sns.gov.pt/api/explore/v2.1/catalog/datasets/evolucao-situacoes-doentes-sinais-sintomas-de-avc/exports/csv?lang=pt&timezone=Europe%2FLondon&use_labels=true&delimiter=%3B")

## Vamos primeiro analisar o dataset hypertension:

#Ver missing values 

out <- sapply(hypertension, function(x) sum(is.na(x) | x == ""))
       
#no dataframe de hipertensão existem 24 missing values na variavel ponto_ou_localizacao_geografica.
#Retiro as observações ou faço imputação? Tenho de ver se é relevante também


#avaliar relevância da variável
#ver relevância da localização para os valores de hipertenção
# Modelo completo - regressão de poisson (variável é uma contagem)
modelo_completo <- glm(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n ~ 
    ponto_ou_localizacao_geografica + regiao + tempo, data = hypertension)

# Modelo sem ponto_ou_localizacao_geografica
modelo_reduzido <- glm(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n ~ 
    regiao + tempo, data = hypertension)

# Teste de razão de verossimilhança
anova(modelo_reduzido, modelo_completo, test = "Chisq")

Analysis of Deviance Table

Model 1: contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n ~ 
    regiao + tempo
Model 2: contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n ~ 
    ponto_ou_localizacao_geografica + regiao + tempo
  Resid. Df Resid. Dev  Df   Deviance  Pr(>Chi)    
1      6432 4.4253e+10                             
2      6317 1.0484e+10 115 3.3769e+10 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#variável relevante


total_data <-  length(hypertension$ponto_ou_localizacao_geografica)
percentage_null <- 100*as.numeric(out["ponto_ou_localizacao_geografica"])/total_data

#percentagem de missing values é 0.3659 - muito pequena logo em princípio pode ser aceitável remove-las


#ver distribuição temporal dos valores sem localização correspondente e respetiva percentagem
where_na <- hypertension %>%
  filter(ponto_ou_localizacao_geografica == "" | is.na (ponto_ou_localizacao_geografica)) %>%
  count(tempo, name = "ausentes")

total_period <- hypertension %>%
  count(tempo, name = "total")

percentages <- left_join(where_na, total_period, by = "tempo") %>%
  mutate(percentagem = round((ausentes / total) * 100, 2))

#a percentagem sobe 5.13% , criando duvidas sobre a remoção

#comportamento dos valores sem localização correspondente
# Adiciona uma flag para missing ou não
hypertension$localizacao <- ifelse(hypertension$ponto_ou_localizacao_geografica == "" | is.na(hypertension$ponto_ou_localizacao_geografica), FALSE, TRUE)

# Compara as médias e medianas
# total de dados
metrica_total <- hypertension %>%
     summarise(across(where(is.numeric), list(media = mean, mediana = median, variancia = var), na.rm = TRUE))
 
# dados com localização
metrica_com_local <- hypertension %>%
     filter(localizacao != FALSE) %>%
     summarise(across(where(is.numeric), list(media = mean, mediana = median, variancia = var), na.rm = TRUE))

#os valores das métricas comparadas não são muito discrepantes logo concluimos que a melhor abordagem é a imputação porque os dados sem localização representam a amostra, a que é confirmado pelo Welch Two Sample t-test
amostra1 <- hypertension$contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n
amostra2 <- hypertension %>%
  filter(localizacao != FALSE) %>%
  pull(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n)

t.test(amostra1, amostra2, var.equal = FALSE)

    Welch Two Sample t-test

data:  amostra1 and amostra2
t = 0.45892, df = 13093, p-value = 0.6463
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -96.10491 154.86345
sample estimates:
mean of x mean of y 
 4945.170  4915.791 

# p-value= 0.6463 > 0.05 , não rejeitamos a hipotese nula. Concluímos que não há diferença significativa entre as médias.

# tornar index da localização mais acessivel para que a imputação possa fazer mais sentido
# Renomear a coluna de coordenadas
hypertension <- hypertension %>%
  rename(coords = ponto_ou_localizacao_geografica)

# Separar coordenadas em lat/lon
# Separar a coluna de coordenadas em latitude e longitude

hyp_coords <- hypertension %>%
  separate(coords, into = c("lat", "lon"), sep = ",", convert = TRUE) %>%
  mutate(
    lat = as.numeric(str_replace(lat, ",", ".")),
    lon = as.numeric(str_replace(lon, ",", ".")),
    proporcao_hipertensos_65_a_com_pa_150_90 = as.numeric(proporcao_hipertensos_65_a_com_pa_150_90)
  )

Warning: Expected 2 pieces. Missing pieces filled with `NA` in 24 rows [1244, 1248, 1253, 1286, 1322,
2503, 2521, 2529, 2537, 2545, 2572, 3880, 3909, 3914, 5207, 5225, 5250, 5262, 5277, 5286, ...].

# Garantir que a proporção esteja em formato numérico
hyp_coords <- hyp_coords %>%
  mutate(proporcao_hipertensos_65_a_com_pa_150_90 = as.numeric(proporcao_hipertensos_65_a_com_pa_150_90))


# Converter para sf
hyp_sf <- hyp_coords %>%
  mutate(geometry = pmap(list(lon, lat), function(x, y) {
    if (is.na(x) || is.na(y)) {
      st_geometrycollection()
    } else {
      st_point(c(x, y))
    }
  })) %>%
  st_as_sf(crs = 4326)
# ----------------------------------
# 4. Carregar shapefiles (níveis 2)
# ----------------------------------

# Municípios (Nível 2)
municipios <- st_read("gadm41_PRT_2.shp") %>% st_transform(4326)

Reading layer `gadm41_PRT_2' from data source 
  `/Users/marianahenriques/Documents/ProjetoBioEst/gadm41_PRT_shp/gadm41_PRT_2.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 308 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -31.26819 ymin: 30.03018 xmax: -6.189142 ymax: 42.15432
Geodetic CRS:  WGS 84

municipios <- municipios %>%
  rename(distritos = NAME_1,
         municipio = NAME_2)

# ----------------------------------
# 5. Join espacial: Hipertensão → Município
# ----------------------------------

hyper_joined <- st_join(hyp_sf, municipios, join = st_within)

#imputacao com base na semelhança de valores

# Calcular a média por grupo municipio
mean_per_group <- hyper_joined %>%
  filter(!is.na(municipio)) %>%
  group_by(municipio) %>%
  summarise(mean_contagem = mean(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n, na.rm = TRUE)) %>%
  st_set_geometry(NULL)

desvio <- hyper_joined %>%
  filter(!is.na(municipio)) %>%
  group_by(municipio) %>%
  summarise(desvio_contagem = sd(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n, na.rm = FALSE)) %>%
  st_set_geometry(NULL)


intervalo <- mean_per_group %>%
  left_join(desvio, by = "municipio") %>%
  mutate(
    limite_inferior = mean_contagem - desvio_contagem,
    limite_superior = mean_contagem + desvio_contagem
  )

# Imputar municipio com base na média mais próxima da contagem
hyp_final <- hyper_joined %>%
  rowwise() %>%
  mutate(municipio = if_else(
    is.na(municipio),
    {
      cont_val <- contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n
      
      # Verifica quais municípios têm o valor dentro do intervalo
      dentro_intervalo <- intervalo %>%
        filter(cont_val >= limite_inferior & cont_val <= limite_superior)

      if (nrow(dentro_intervalo) > 0) {
        # Entre os válidos, pega o de média mais próxima
        diffs <- abs(dentro_intervalo$mean_contagem - cont_val)
        dentro_intervalo$municipio[which.min(diffs)]
      } else {
        NA_character_
      }
    },
    municipio
  )) %>%
  ungroup()%>%
  filter(!is.na(municipio))



tabela_correspondencia <- municipios %>%
  st_drop_geometry() %>%
  select(distritos, municipio) %>%
  distinct()

hyp_final <- hyp_final %>%
  mutate(municipio_clean = stri_trans_general(tolower(trimws(municipio)), "Latin-ASCII"))

tabela_correspondencia <- tabela_correspondencia %>%
  mutate(municipio_clean = stri_trans_general(tolower(trimws(municipio)), "Latin-ASCII"))

# 2. Faz o join com base na versão limpa
hyp_final <- hyp_final %>%
  left_join(tabela_correspondencia %>% select(municipio_clean, distritos), by = "municipio_clean")

# --- Efeito da imputação na média ---
# 1. Média original por grupo
mean_attr <- hyper_joined %>%
  filter(!is.na(municipio)) %>%
  group_by(municipio) %>%
  summarise(mean_original = mean(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n, na.rm = TRUE)) %>%
  st_set_geometry(NULL)

# 2. Média após imputação (em hyp_final)
mean_per_group_com_imp <- hyp_final %>%
  group_by(municipio) %>%
  summarise(mean_imputada = mean(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n, na.rm = TRUE)) %>%
  st_set_geometry(NULL)

# 3. Juntar médias original e imputada
means <- left_join(mean_attr, mean_per_group_com_imp, by = "municipio")

# 4. Gráfico das médias
means_long <- means %>%
  pivot_longer(cols = c(mean_original, mean_imputada),
               names_to = "tipo_media",
               values_to = "media")

grafico <- ggplot(means_long, aes(x = municipio, y = media, fill = tipo_media)) +
  geom_col(position = "dodge") +
  labs(title = "Comparação das médias por municipio",
       x = "Municipio", y = "Média da contagem",
       fill = "Tipo de Média") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggsave(filename = "grafico_medias_por_grupo.png", plot = grafico, width = 10, height = 6, dpi = 300)


# 5. Distância entre valor imputado e média original do grupo
distance <- hyp_final %>%
  filter(localizacao == FALSE) %>%
  st_set_geometry(NULL) %>%
  left_join(mean_attr, by = "municipio") %>%
  mutate(dist = mean_original - contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n)

# 6. Gráfico de linhas para analisar 5.
# Criar limites para a região sombreada (banda de confiança)
distance <- distance %>%
  mutate(
    ymin = mean_original - abs(dist),
    ymax = mean_original + abs(dist)
  )

grafico_2 <- ggplot(distance, aes(x = reorder(municipio, mean_original), y = mean_original, group = 1)) +
  geom_ribbon(aes(ymin = ymin, ymax = ymax), fill = "lightblue", alpha = 0.4) +  # região sombreada
  geom_line(color = "blue", size = 1) +  # linha da média
  geom_point(aes(y = contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n),
             color = "red", size = 2) +  # pontos do valor imputado
  labs(
    x = "Municípios com valores imputados",
    y = "Valor",
    title = "Média original vs distância do valor imputado",
    subtitle = "Linha azul: média original; Região azul clara: distância; Pontos vermelhos: valores imputados"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

ggsave(filename = "grafico_distancias_por_grupo.png", plot = grafico_2, width = 10, height = 6, dpi = 300)

# --- Comparação da distribuição regional: imputados vs reais ---
tabela_comparativa <- hyp_final %>%
  mutate(tipo = ifelse(localizacao, "Real", "Imputado")) %>%
  group_by(municipio, tipo) %>%
  summarise(total = n(), .groups = "drop") %>%
  group_by(tipo) %>%
  mutate(prop = total / sum(total))

p <- ggplot(tabela_comparativa, aes(x = prop, y = municipio, fill = tipo)) +
  geom_col(position = "dodge") +
  theme_minimal() +
  labs(title = "Distribuição de Municípios: Real vs Imputado",
       y = "Município", x = "Proporção") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

ggsave("grafico_distribuicao.png", plot = p, width = 10, height = 8, dpi = 300)

# --- Correlação ---
cor_before <- cor(
  hyper_joined$contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n,
  hyper_joined$proporcao_hipertensos_65_a_com_pa_150_90,
  use = "complete.obs"
)

cor_after <- cor(
  hyp_final$contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n,
  hyp_final$proporcao_hipertensos_65_a_com_pa_150_90
)

cat("Correlação antes da imputação:", cor_before, "\n")

Correlação antes da imputação: 0.7432419 

cat("Correlação depois da imputação:", cor_after, "\n")

Correlação depois da imputação: 0.7427319 

#summary
summary(hyper_joined %>% select(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n, proporcao_hipertensos_65_a_com_pa_150_90))

 contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n
 Min.   :  129                                                                                       
 1st Qu.: 2273                                                                                       
 Median : 4117                                                                                       
 Mean   : 4945                                                                                       
 3rd Qu.: 6630                                                                                       
 Max.   :43441                                                                                       
 proporcao_hipertensos_65_a_com_pa_150_90               geometry   
 Min.   : 1.32                            GEOMETRYCOLLECTION:  24  
 1st Qu.:18.07                            POINT             :6536  
 Median :29.83                            epsg:4326         :   0  
 Mean   :31.99                            +proj=long...     :   0  
 3rd Qu.:43.24                                                     
 Max.   :81.75                                                     

summary(hyp_final %>% select(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n, proporcao_hipertensos_65_a_com_pa_150_90))

 contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n
 Min.   :  129                                                                                       
 1st Qu.: 2272                                                                                       
 Median : 4115                                                                                       
 Mean   : 4935                                                                                       
 3rd Qu.: 6627                                                                                       
 Max.   :43441                                                                                       
 proporcao_hipertensos_65_a_com_pa_150_90               geometry   
 Min.   : 1.32                            GEOMETRYCOLLECTION:  20  
 1st Qu.:18.07                            POINT             :6536  
 Median :29.82                            epsg:4326         :   0  
 Mean   :31.97                            +proj=long...     :   0  
 3rd Qu.:43.18                                                     
 Max.   :81.75                                                     

# --- Dispersão das variáveis principais ---
ggplot(hyp_final, aes(x = contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n,
                      y = proporcao_hipertensos_65_a_com_pa_150_90,
                      color = localizacao)) +
  geom_point() +
  labs(title = "Correlação: Reais vs Imputados", color = "Dados reais?")

# --- PCA para visualização multivariada ---
dados_pca <- st_drop_geometry(hyp_final[, c("contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n","proporcao_hipertensos_65_a_com_pa_150_90")])

res.pca <- PCA(dados_pca, graph = FALSE)
fviz_pca_ind(res.pca, 
             habillage = as.factor(hyp_final$localizacao),
             title = "PCA: Imputed vs Real Data",
             palette = c("red", "blue"))  # Vermelho = imputado, Azul = real

ggplot(hyper_joined, aes(x = contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.5) + ggtitle("Contagem Antes da Imputação")

ggplot(hyp_final, aes(x = contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n)) +
  geom_histogram(bins = 30, fill = "green", alpha = 0.5) + ggtitle("Contagem Depois da Imputação")

ggplot(hyper_joined, aes(x = proporcao_hipertensos_65_a_com_pa_150_90)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.5) + ggtitle("Proporção Antes da Imputação")

ggplot(hyp_final, aes(x = proporcao_hipertensos_65_a_com_pa_150_90)) +
  geom_histogram(bins = 30, fill = "green", alpha = 0.5) + ggtitle("Proporção Depois da Imputação")

#imputação com sucesso

#Analise temporal do dataframe da hipertensao com imputacao

# Calcular proporção total
hyp_final <- hyp_final %>%
  mutate(total = 100 * contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n / proporcao_hipertensos_65_a_com_pa_150_90)

# Total por ano (todos os utentes)
hipertensao_ano <- hyp_final %>%
  mutate(Ano = year(ym(tempo))) %>%
  group_by(Ano) %>%
  summarise(total = sum(as.numeric(total), na.rm = TRUE), .groups = "drop") %>%
  st_set_geometry(NULL)

ggplot(hipertensao_ano, aes(x = factor(Ano), y = total)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(total)), vjust = -0.25, size = 3) +
  labs(title = "Total de Utentes com Hipertensão por Ano",
       x = "Ano",
       y = "Total de Utentes") +
  theme_minimal()

# Total por ano para pessoas com menos de 65 anos
hipertensao_ano_menos_65 <- hyp_final %>%
  mutate(Ano = year(ym(tempo))) %>%
  group_by(Ano) %>%
  summarise(hip_menos_65 = sum(as.numeric(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n), na.rm = TRUE), .groups = "drop") %>%
  st_set_geometry(NULL)

ggplot(hipertensao_ano_menos_65, aes(x = factor(Ano), y = hip_menos_65)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(hip_menos_65)), vjust = -0.25, size = 3) +
  labs(title = "Total de Utentes com menos de 65 anos com Hipertensão por Ano",
       x = "Ano",
       y = "Total de Utentes") +
  theme_minimal()

# Comparação entre menos de 65 e 65 ou mais
hipertensao_comp_ano <- hipertensao_ano %>%
  left_join(hipertensao_ano_menos_65, by = "Ano") %>%
  mutate(
    hip_menos_65 = ifelse(is.na(hip_menos_65), 0, hip_menos_65),
    hip_mais_65 = total - hip_menos_65
  ) %>%
  select(Ano, hip_menos_65, hip_mais_65) %>%
  pivot_longer(
    cols = c(hip_menos_65, hip_mais_65),
    names_to = "Group",
    values_to = "total"
  ) %>%
  mutate(
    Group = dplyr::recode(Group,
                          "hip_menos_65" = "Under 65 years old",
                          "hip_mais_65" = "65 years or older")
  )

# Gráfico empilhado
ggplot(hipertensao_comp_ano, aes(x = factor(Ano), y = total, fill = Group)) +
  geom_col() +
  labs(
    title = " Controlled Hypertension by Age and Year",
    x = "Year",
    y = "Number of Patients"
  ) +
  theme_minimal()

# Hipertensão por mês - Total
hipertensao_mes <- hyp_final %>%
  mutate(Data = ym(tempo),
         Mes = month(Data, label = TRUE, abbr = FALSE)) %>%  
  group_by(Mes) %>%
  summarise(total = sum(as.numeric(total), na.rm = TRUE), .groups = "drop") %>%
  arrange(match(Mes, month.name)) %>%
  st_set_geometry(NULL)

ggplot(hipertensao_mes, aes(x = Mes, y = total)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(total)), vjust = -0.25, size = 3) +
  labs(title = "Total de Utentes com Hipertensão por Mês",
       x = "Mês",
       y = "Total de Utentes") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Hipertensão por mês - Menos de 65
hipertensao_mes_65 <- hyp_final %>%
  mutate(Data_65 = ym(tempo),
         Mes_65 = month(Data_65, label = TRUE, abbr = FALSE)) %>%  
  group_by(Mes_65) %>%
  summarise(hip_menos_65 = sum(as.numeric(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n), na.rm = TRUE), .groups = "drop") %>%
  arrange(match(Mes_65, month.name)) %>%
  st_set_geometry(NULL)

ggplot(hipertensao_mes_65, aes(x = Mes_65, y = hip_menos_65)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(hip_menos_65)), vjust = -0.25, size = 3) +
  labs(title = "Utentes <65 anos com Hipertensão por Mês",
       x = "Mês",
       y = "Total de Utentes") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Tabela de meses em português
meses_pt <- c("janeiro", "fevereiro", "março", "abril", "maio", "junho", 
              "julho", "agosto", "setembro", "outubro", "novembro", "dezembro")


mes_lookup <- c(
  "january" = "janeiro", "february" = "fevereiro", "march" = "março",
  "april" = "abril", "may" = "maio", "june" = "junho",
  "july" = "julho", "august" = "agosto", "september" = "setembro",
  "october" = "outubro", "november" = "novembro", "december" = "dezembro"
)


# Juntar total e <65, calcular >=65
hipertensao_comp_mes <- hipertensao_mes %>%
  left_join(hipertensao_mes_65, by = c("Mes" = "Mes_65")) %>%
  mutate(
    hip_menos_65 = ifelse(is.na(hip_menos_65), 0, hip_menos_65),
    hip_mais_65 = total - hip_menos_65,
    Mes = dplyr::recode(tolower(as.character(Mes)), !!!mes_lookup),
    Mes = factor(Mes, levels = meses_pt)
  ) %>%
  select(Mes, hip_menos_65, hip_mais_65) %>%
  pivot_longer(
    cols = c(hip_menos_65, hip_mais_65),
    names_to = "Group",
    values_to = "Total_hip"
  ) %>%
  mutate(
    Group = dplyr::recode(Group,
                          "hip_menos_65" = "Under 65 years old",
                          "hip_mais_65" = "65 years or older")
  )

# Gráfico final - barras egroup_by()# Gráfico final - barras empilhadas
ggplot(hipertensao_comp_mes, aes(x = Mes, y = Total_hip, fill = Group)) +
  geom_col() +
  labs(
    title = " Controlled Hypertension by Age and Month",
    x = "Month",
    y = "Number of patients"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

## Análise do dataset avc: 

#Quais os anos e meses é que estão registados no AVC
meses_por_ano_avc <- avc %>%
  mutate(Ano_avc = year(ym(Período)),
         Mes_avc = month(ym(Período), label = TRUE, abbr = FALSE)) %>%
  distinct(Ano_avc, Mes_avc) %>%
  arrange(Ano_avc, Mes_avc) %>%
  group_by(Ano_avc) %>%
  summarise(Meses_Disponiveis = paste(Mes_avc, collapse = ", "))

#Ver missing values 
sapply(avc, function(x) sum(is.na(x) | x == ""))

                 N.º.Registos                       Período        Distrito.da.Ocorrência 
                            0                             0                             0 
       Concelho.da.Ocorrência        Faixa.Etária.da.Vítima              Género.da.Vítima 
                            0                             0                             0 
          Hospital.de.Destino N.º.de.registos.Via.Verde.AVC 
                            0                             0 

# Distribuição por Género com AVC 

avc_genero <- avc %>%
  group_by(Género.da.Vítima) %>%
  summarise(Total_AVC = sum(as.numeric(`N.º.de.registos.Via.Verde.AVC`), na.rm = TRUE)) %>%
  arrange(desc(Total_AVC))

ggplot(avc_genero, aes(x = reorder(Género.da.Vítima, -Total_AVC), y = Total_AVC, fill = Género.da.Vítima)) +
  geom_col() +
  geom_text(aes(label = Total_AVC), vjust = -0.2, size = 2) +
  labs(title = "Total de Registos de AVC por Género",
       x = "Género",
       y = "Total de AVC") +
  theme_minimal()

# Distribuição por Género com AVC em cada ano

avc_genero_ano <- avc %>%
  mutate(Ano = year(ym(Período))) %>%
  group_by(Ano, Género.da.Vítima) %>%
  summarise(Total_AVC = sum(as.numeric(N.º.de.registos.Via.Verde.AVC), na.rm = TRUE), .groups = "drop")

# Traduzir valores de 'Género.da.Vítima'
avc$Género.da.Vítima <- ifelse(avc$Género.da.Vítima == "Feminino", "Female",
                               ifelse(avc$Género.da.Vítima == "Masculino", "Male", avc$Género.da.Vítima))

ggplot(avc_genero_ano, aes(x = factor(Ano), y = Total_AVC, fill = Género.da.Vítima)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = Total_AVC), position = position_dodge(width = 0.9), vjust = -0.25, size = 3) +
  labs(title = "Total Stroke Records by Year and Gender",
       x = "Year",
       y = "Total Stroke Records",
       fill = "Gender") +
  theme_minimal()

# Distribuição por Género com AVC em cada mês
avc_genero_mes <- avc %>%
  mutate(Mes_Ano = format(ym(Período), "%Y-%m")) %>%
  group_by(Mes_Ano, `Género.da.Vítima`) %>%
  summarise(Total_AVC = sum(as.numeric(`N.º.de.registos.Via.Verde.AVC`), na.rm = TRUE), .groups = "drop")

ggplot(avc_genero_mes, aes(x = Mes_Ano, y = Total_AVC, fill = `Género.da.Vítima`)) +
  geom_col(position = "dodge") +
  labs(title = "Stroke Distribution by Month and Gender",
       x = "Month", y = "Total Stroke Records", fill = "Gender") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Distribuição Faixa Etária com AVC no total dos anos

# Define the correct order of age ranges
faixa_ordem <- c("0 - 17", "18 - 29", "30 - 39", "40 - 49", "50 - 59", 
                 "60 - 69", "70 - 79", "80 - 89", "90 - 99", ">100")

# Apply the order using factor
avc_idade <- avc %>%
  group_by(`Faixa.Etária.da.Vítima`) %>%
  summarise(Total_AVC = sum(as.numeric(`N.º.de.registos.Via.Verde.AVC`), na.rm = TRUE)) %>%
  mutate(`Faixa.Etária.da.Vítima` = factor(`Faixa.Etária.da.Vítima`, levels = faixa_ordem)) %>%
  arrange(`Faixa.Etária.da.Vítima`)

# Plot
ggplot(avc_idade, aes(x = `Faixa.Etária.da.Vítima`, y = Total_AVC)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = Total_AVC), vjust = -0.2, size = 2) +
  labs(title = "Total de Registos de AVC por Faixa Etária",
       x = "Faixa Etária",
       y = "Total de AVC") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Distribuição Faixa Etária com AVC em cada ano

avc_faixa_ano <- avc %>%
  mutate(Ano = year(ym(Período))) %>%
  group_by(Ano, `Faixa.Etária.da.Vítima`) %>%
  summarise(Total_AVC = sum(as.numeric(`N.º.de.registos.Via.Verde.AVC`), na.rm = TRUE), .groups = "drop") %>% mutate(`Faixa.Etária.da.Vítima` = factor(`Faixa.Etária.da.Vítima`, levels = faixa_ordem)) %>%
  arrange(`Faixa.Etária.da.Vítima`)


ggplot(avc_faixa_ano, aes(x = factor(Ano), y = Total_AVC, fill = `Faixa.Etária.da.Vítima`)) +
  geom_col(position = "dodge") +
  labs(title = "Total Stroke Records by Year and Age Group",
       x = "Year",
       y = "Total Stroke Records",
       fill = "Age Group") +
  theme_minimal()

# Distribuição Faixa Etária com AVC em cada mês

avc_faixa_mes <- avc %>%
  mutate(Mes_Ano = format(ym(Período), "%Y-%m")) %>%
  group_by(Mes_Ano, `Faixa.Etária.da.Vítima`) %>%
  summarise(Total_AVC = sum(as.numeric(`N.º.de.registos.Via.Verde.AVC`), na.rm = TRUE), .groups = "drop") %>% mutate(`Faixa.Etária.da.Vítima` = factor(`Faixa.Etária.da.Vítima`, levels = faixa_ordem)) %>%
  arrange(`Faixa.Etária.da.Vítima`)

ggplot(avc_faixa_mes, aes(x = Mes_Ano, y = Total_AVC, fill = `Faixa.Etária.da.Vítima`)) +
  geom_col(position = "dodge") +
  labs(title = "Stroke Distribution by Month and Age Group",
       x = "Month", y = "Total Stroke Records", fill = "Age Groupe") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Género × Faixa Etária ao longos do período 

avc_cross <- avc %>%
  group_by(`Género.da.Vítima`, `Faixa.Etária.da.Vítima`) %>%
  summarise(Total_AVC = sum(as.numeric(`N.º.de.registos.Via.Verde.AVC`), na.rm = TRUE)) %>%
  ungroup() %>% mutate(`Faixa.Etária.da.Vítima` = factor(`Faixa.Etária.da.Vítima`, levels = faixa_ordem)) %>%
  arrange(`Faixa.Etária.da.Vítima`)

`summarise()` has grouped output by 'Género.da.Vítima'. You can override using the `.groups`
argument.

ggplot(avc_cross, aes(x = `Faixa.Etária.da.Vítima`, y = Total_AVC, fill = `Género.da.Vítima`)) +
  geom_col(position = "dodge") +
  labs(title = "Distribuição de AVC por Género e Faixa Etária", x = "Faixa Etária", y = "Total de AVC", fill = "Género") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Correlação entre Hipertensos e AVC por grupo etário (maiores de 65 anos e menores de 65 anos)

# Clean names
avc <- avc %>% janitor::clean_names()

# Fix column name
names(avc)[names(avc) == "n_o_de_registos_via_verde_avc"] <- "n_de_registos_via_verde_avc"

# Dividimos o AVC em maiores de 65 anos e menores de 65 anos
avc_clean <- avc %>%
  mutate(
    faixa_inicio = as.numeric(str_extract(faixa_etaria_da_vitima, "^\\d+")),
    faixa_fim = as.numeric(str_extract(faixa_etaria_da_vitima, "\\d+$")),
    Group = case_when(
      faixa_inicio < 69 ~ "Under 65 years old",
      TRUE ~ "65 years or older"
    ),
    mes_num = as.integer(str_sub(periodo, 6, 7)),
    mes = factor(meses_pt[mes_num], levels = meses_pt),
    n_de_registos_via_verde_avc = as.numeric(n_de_registos_via_verde_avc)
  ) %>%
  group_by(mes, Group) %>%
  summarise(total_avc = sum(n_de_registos_via_verde_avc, na.rm = TRUE), .groups = "drop")

# Filtrar dados de hipertensão de 2022-01 a 2024-02
hyp_filtrado <- hyp_final %>%
  st_set_geometry(NULL) %>%
  mutate(Data = ym(tempo)) %>%
  filter(Data >= ym("2022-01") & Data <= ym("2024-02"))

# Hipertensão total por mês
hipertensao_mes_filtrado <- hyp_filtrado %>%
  mutate(Mes = month(Data, label = TRUE, abbr = FALSE)) %>%
  group_by(Mes) %>%
  summarise(total = sum(as.numeric(total), na.rm = TRUE), .groups = "drop") %>%
  arrange(match(Mes, month.name))

# Hipertensão em menores de 65 anos
hipertensao_mes_65_filtrado <- hyp_filtrado %>%
  mutate(Mes = month(Data, label = TRUE, abbr = FALSE)) %>%
  group_by(Mes) %>%
  summarise(hip_65 = sum(as.numeric(contagem_de_utentes_inscritos_com_hipertensao_arterial_com_pressao_arterial_inferior_a_150_90_mmhg_n), na.rm = TRUE), .groups = "drop") %>%
  arrange(match(Mes, month.name))

# Combinar dados de hipertensão e preparar formato longo
hipertensao_comp_mes_filtrado <- hipertensao_mes_filtrado %>%
  left_join(hipertensao_mes_65_filtrado, by = "Mes") %>%
  mutate(
    hip_65 = ifelse(is.na(hip_65), 0, hip_65),
    hip_mais_65 = total - hip_65,
    Mes = mes_lookup[Mes],
    Mes = factor(Mes, levels = meses_pt)
  ) %>%
  select(Mes, hip_65, hip_mais_65) %>%
  pivot_longer(
    cols = c("hip_65", "hip_mais_65"),
    names_to = "Group",
    values_to = "Total_hip"
  ) %>%
  mutate(
    Group = dplyr::recode(Group,
      "hip_65" = "Under 65 years old",
      "hip_mais_65" = "65 years or older")
  )

# Gráfico empilhado de hipertensão
ggplot(hipertensao_comp_mes_filtrado, aes(x = Mes, y = Total_hip, fill = Group)) +
  geom_col() +
  labs(
    title = "Hipertensão por Idade e Mês (2022-01 a 2024-02)",
    x = "Mês",
    y = "Número de Utentes"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Resumo para análise
hipertensao_summary <- hipertensao_comp_mes_filtrado %>%
  group_by(Mes, Group) %>%
  summarise(Total_HIP = sum(Total_hip), .groups = "drop") %>%
  rename(mes = Mes)

# Merge com dados de AVC
merged <- left_join(hipertensao_summary, avc_clean, by = c("mes", "Group"))

# Boxplots para visualização de outliers
ggplot(merged, aes(x = Group, y = Total_HIP)) +
  geom_boxplot() +
  labs(title = "Boxplot of Controlled Hypertension by Age Group", x="Group", y="Total Hypertension Patients")

ggplot(merged, aes(x = Group, y = total_avc)) +
  geom_boxplot() +
  labs(title = "Boxplot of Stroke Cases by Age Group", x="Group", y="Total Stroke Cases")

#como há outliers, fazemos com o spearman

# Correlação de Spearman
cor_results_spearman <- merged %>%
  group_by(Group) %>%
  summarise(correl = cor(Total_HIP, total_avc, method = "spearman"))

print(cor_results_spearman)

# Scatterplot com regressão para visualização
ggplot(merged, aes(x = Total_HIP, y = total_avc, color = Group)) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ Group) +
  labs(
    title = "Correlation between Controlled Hypertension and Stroke by Age Group",
    x = "Total Hypertension Patients",
    y = "Total Stroke Cases"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

# Regressões lineares por grupo
models <- merged %>%
  group_by(Group) %>%
  group_map(~ lm(total_avc ~ Total_HIP, data = .x), .keep = TRUE)

# Resumo das regressões
summary(models[[1]])  # Under 65 years old

Call:
lm(formula = total_avc ~ Total_HIP, data = .x)

Residuals:
    Min      1Q  Median      3Q     Max 
-152.51  -36.62   16.85   56.54  106.31 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.703e+02  7.398e+01   6.357 8.28e-05 ***
Total_HIP   3.116e-04  5.918e-05   5.266 0.000365 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 82.79 on 10 degrees of freedom
Multiple R-squared:  0.735, Adjusted R-squared:  0.7085 
F-statistic: 27.73 on 1 and 10 DF,  p-value: 0.0003649

summary(models[[2]])  # 65 years or older

Call:
lm(formula = total_avc ~ Total_HIP, data = .x)

Residuals:
   Min     1Q Median     3Q    Max 
-83.16 -29.66 -13.39  10.95 148.96 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.935e+02  5.092e+01   9.691 2.12e-06 ***
Total_HIP   -1.397e-04  8.325e-05  -1.679    0.124    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 64.71 on 10 degrees of freedom
Multiple R-squared:  0.2198,    Adjusted R-squared:  0.1418 
F-statistic: 2.817 on 1 and 10 DF,  p-value: 0.1242

# Total de Registos de AVC por Distrito
avc_regiao <- avc %>%
  group_by(`distrito_da_ocorrencia`) %>%
  summarise(Total_AVC = sum(as.numeric(`n_de_registos_via_verde_avc`), na.rm = TRUE)) %>%
  arrange(desc(Total_AVC)) 

ggplot(avc_regiao, aes(x = reorder(`distrito_da_ocorrencia`, -Total_AVC), y = Total_AVC)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = Total_AVC), vjust = -0.2, size = 2) +
  labs(title = "Total Stroke Records by District",
       x = "District", y = "Patients") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Distribuição ao longo do tempo por distrito

avc_regiao_tempo <- avc %>%
  mutate(Mes_Ano = format(ym(periodo), "%Y-%m")) %>%
  group_by(Mes_Ano, `distrito_da_ocorrencia`) %>%
  summarise(Total_AVC = sum(as.numeric(`n_de_registos_via_verde_avc`), na.rm = TRUE), .groups = "drop")

ggplot(avc_regiao_tempo, aes(x = Mes_Ano, y = Total_AVC, color = `distrito_da_ocorrencia`, group = `distrito_da_ocorrencia`)) +
  geom_line() +
  labs(title = "Evolução Mensal de Registos de AVC por Distrito",
       x = "Mês", y = "Total de AVC", color = "Distrito") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## Ir buscar o número de habitantes por distrito em portugal em 2022, 2023 e 2024 e fazer média, para normalizar e poder comparar, porque obviamente que no porto e lisboa há maior numero de avc porque há mais habitantes.

# Vetor com as populações médias estimadas
populacao_distritos <- c(
  "Lisboa" = 2318952,
  "Porto" = 1829758,
  "Setúbal" = 893464,
  "Braga" = 858708,
  "Aveiro" = 717998,
  "Faro" = 474500,
  "Leiria" = 471570,
  "Santarém" = 436891,
  "Coimbra" = 413225,
  "Viseu" = 364020,
  "Madeira" = 254630,
  "Açores" = 241471,
  "Viana do Castelo" = 233610,
  "Vila Real" = 185263,
  "Castelo Branco" = 178860,
  "Évora" = 153427,
  "Beja" = 147363,
  "Guarda" = 142131,
  "Bragança" = 123109,
  "Portalegre" = 104561
)


# Adicionar população ao DataFrame
avc_regiao_tempo <- avc_regiao_tempo %>%
  mutate(Populacao = populacao_distritos[`distrito_da_ocorrencia`],
         Taxa_AVC = (Total_AVC / Populacao) * 100000)

# Média da taxa de AVC por distrito
ranking_distritos <- avc_regiao_tempo %>%
  group_by(`distrito_da_ocorrencia`) %>%
  summarise(Media_Taxa_AVC = mean(Taxa_AVC, na.rm = TRUE)) %>%
  arrange(desc(Media_Taxa_AVC))

print(ranking_distritos)

ggplot(ranking_distritos, aes(x = reorder(`distrito_da_ocorrencia`, Media_Taxa_AVC), y = Media_Taxa_AVC)) +
  geom_col(fill = "darkred") +
  geom_text(aes(label = round(Media_Taxa_AVC, 1)), hjust = -0.1, size = 3) +
  coord_flip() +
  labs(title = "Stroke Rate by District",
       x = "District", y = "Rate") +
  theme_minimal()

# Visualizar as taxas de AVC ao longo do tempo por distrito
ggplot(avc_regiao_tempo, aes(x = Mes_Ano, y = Taxa_AVC, color = `distrito_da_ocorrencia`, group = `distrito_da_ocorrencia`)) +
  geom_line() +
  labs(title = "Evolução Mensal da Taxa de AVC por Distrito",
       x = "Mês", y = "Taxa de AVC por 100.000 habitantes", color = "Distrito") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

#Como não é possivel comparar bem pelo gráfico, vamos fazer testes estatitisticos. 

#vamos ver primeiro se o total_avc tem distribuição normal em todos os distritos

ggplot(avc_regiao_tempo, aes(x = Taxa_AVC)) +
  geom_histogram(bins = 15, fill = "skyblue", color = "black") +
  facet_wrap(~ `distrito_da_ocorrencia`, scales = "free") +
  theme_minimal()

# Aplicar o teste shapiro 
shapiro_results <- avc_regiao_tempo %>%
  group_by(distrito_da_ocorrencia) %>%
  summarise(p_value = shapiro.test(Taxa_AVC)$p.value)

# Print
print(shapiro_results)

# Homogeneidade das variâncias entre os grupos

levene_result <- leveneTest(Taxa_AVC ~ distrito_da_ocorrencia, data = avc_regiao_tempo)

Warning in leveneTest.default(y = y, group = group, ...) :
  group coerced to factor.

# Printar o resultado
print(levene_result)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value    Pr(>F)    
group  17  9.8049 < 2.2e-16 ***
      450                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#conclusão: os dados são distribuidos normalmente e as variâncias são significativamente diferentes. Logo, vamos fazer o teste Welch ANOVA (oneway.test) para comparar as médias das taxas de AVC por 100.000 habitantes entre distritos, sem assumir variâncias iguais.

# Welch ANOVA
oneway.test(Taxa_AVC ~ `distrito_da_ocorrencia`, data = avc_regiao_tempo, var.equal = FALSE)

    One-way analysis of means (not assuming equal variances)

data:  Taxa_AVC and distrito_da_ocorrencia
F = 32.86, num df = 17.00, denom df = 167.03, p-value < 2.2e-16

#vamos ver como a taxa média de AVC evoluiu ao longo dos anos por distrito

avc_regiao_tempo <- avc_regiao_tempo %>%
  mutate(Ano = substr(Mes_Ano, 1, 4))  # extrai os 4 primeiros caracteres como ano

media_anual_distrito <- avc_regiao_tempo %>%
  group_by(Ano, `distrito_da_ocorrencia`) %>%
  summarise(Media_Taxa_AVC = mean(Taxa_AVC, na.rm = TRUE), .groups = "drop")

ggplot(media_anual_distrito, aes(x = Ano, y = Media_Taxa_AVC, color = `distrito_da_ocorrencia`, group = `distrito_da_ocorrencia`)) +
  geom_line(size = 1) +
  geom_point(size = 1.5) +
  labs(title = "Evolução Anual da Taxa de AVC por Distrito",
       x = "Ano",
       y = "Taxa média de AVC por 100.000 habitantes",
       color = "Distrito") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

#  Heatmap 
ggplot(media_anual_distrito, aes(x = Ano, y = reorder(`distrito_da_ocorrencia`, -Media_Taxa_AVC))) +
  geom_tile(aes(fill = Media_Taxa_AVC), color = "white") +
  scale_fill_gradient(low = "white", high = "red") +
  labs(title = "Heat Map: Stroke Rate by District and Year",
       x = "Year",
       y = "District",
       fill = "Rate") +
  theme_minimal()

## A mesma análise mas para hipertensos 
# 1. Dados mensais filtrados (2022-01 a 2024-02)
hipertensao_tempo <- hyp_final %>%
  mutate(Mes_Ano = format(ym(tempo), "%Y-%m")) %>%
  filter(Mes_Ano >= "2022-01", Mes_Ano <= "2024-02") %>%
  group_by(Mes_Ano, distritos.x) %>%
  summarise(
    Total_Hipertensao_Controlada = sum(as.numeric(total),na.rm = TRUE),
    .groups = "drop") %>%
  mutate( Populacao = populacao_distritos[distritos.x],
          Taxa_Hipertensao_Controlada = (Total_Hipertensao_Controlada / Populacao) * 100000
  )%>%
  st_drop_geometry(hipertensao_tempo)%>%
  rename(Distrito = distritos.x)

# 2. Total acumulado por distrito
hipertensao_total <- hipertensao_tempo %>%
  group_by(Distrito) %>%
  summarise(Total_Hipertensao_Controlada = sum(Total_Hipertensao_Controlada, na.rm = TRUE)) %>%
  arrange(desc(Total_Hipertensao_Controlada))


# Gráfico de total acumulado
ggplot(hipertensao_total, aes(x = reorder(Distrito, -Total_Hipertensao_Controlada), y = Total_Hipertensao_Controlada)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = round(Total_Hipertensao_Controlada)), vjust = -0.2, size = 2) +
  labs(title = "Total Patients with Controlled Hypertension by District",
       x = "District", y = "Patients") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Distribuição ao longo do tempo por distrito

ggplot(hipertensao_tempo, aes(x = Mes_Ano, y = Total_Hipertensao_Controlada, color = `Distrito`, group = `Distrito`)) +
  geom_line() +
  labs(title = "Evolução Mensal de Registos de hipertensao por Distrito",
       x = "Mês", y = "Total de Hipertensão", color = "Distrito") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 3. Gráfico de evolução mensal das taxas
ggplot(hipertensao_tempo, aes(x = Mes_Ano, y = Taxa_Hipertensao_Controlada, color = Distrito, group = Distrito)) +
  geom_line() +
  labs(title = "Evolução mensal da Taxa de Hipertensão Controlada por Distrito",
       x = "Mês", y = "Taxa de Hipertensão por 100.000 habitantes", color = "Distrito") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1
  ))

# 4. Média da taxa por distrito
ranking_distritos <- hipertensao_tempo %>%
  group_by(Distrito) %>%
  summarise(Media_Taxa_Hipertensao_Controlada = mean(Taxa_Hipertensao_Controlada, na.rm = TRUE)) %>%
  arrange(desc(Media_Taxa_Hipertensao_Controlada))

# Gráfico de ranking das taxas médias
ggplot(ranking_distritos, aes(x = reorder(Distrito, Media_Taxa_Hipertensao_Controlada), y = Media_Taxa_Hipertensao_Controlada)) +
  geom_col(fill = "darkgreen") +
  geom_text(aes(label = round(Media_Taxa_Hipertensao_Controlada, 1)), hjust = -0.1, size = 2) +
  coord_flip() +
  labs(title = "Controlled Hypertension Rate by District",
       x = "District", y = "Rate") +
  theme_minimal()

# 5. Histogramas por distrito (distribuição das taxas)
ggplot(hipertensao_tempo, aes(x = Taxa_Hipertensao_Controlada)) +
  geom_histogram(bins = 15, fill = "skyblue", color = "black") +
  facet_wrap(~ Distrito, scales = "free") +
  theme_minimal()

# 6. Teste de normalidade (Shapiro-Wilk)
shapiro_results <- hipertensao_tempo %>%
  group_by(Distrito) %>%
  summarise(p_value = shapiro.test(Taxa_Hipertensao_Controlada)$p.value)

print(shapiro_results)

# 7. Teste de homogeneidade das variâncias (Levene)
leveneTest(Taxa_Hipertensao_Controlada ~ Distrito, data = hipertensao_tempo)

Warning in leveneTest.default(y = y, group = group, ...) :
  group coerced to factor.

Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)  
group  17  1.9335 0.01402 *
      449                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 
kruskal.test(Taxa_Hipertensao_Controlada ~ Distrito, data = hipertensao_tempo)

    Kruskal-Wallis rank sum test

data:  Taxa_Hipertensao_Controlada by Distrito
Kruskal-Wallis chi-squared = 406.15, df = 17, p-value < 2.2e-16

# Partindo do objeto `hipertensao_tempo` já existente e filtrado (2022-2024)

# 1. Extrair o ano
hipertensao_tempo <- hipertensao_tempo %>%
  mutate(Ano = substr(Mes_Ano, 1, 4))

# 2. Calcular média anual da taxa por distrito
media_anual_distrito <- hipertensao_tempo %>%
  group_by(Ano, Distrito) %>%
  summarise(Media_Taxa_Hipertensao_Controlada = mean(Taxa_Hipertensao_Controlada, na.rm = TRUE), .groups = "drop")

# 3. Gráfico de linhas: Evolução anual
ggplot(media_anual_distrito, aes(x = Ano, y = Media_Taxa_Hipertensao_Controlada, color = Distrito, group = Distrito)) +
  geom_line(size = 1) +
  geom_point(size = 1.5) +
  labs(title = "Evolução Anual da Taxa de Hipertensão Controlada por Distrito",
       x = "Ano",
       y = "Taxa média por 100.000 habitantes",
       color = "Distrito") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 4. Heatmap: Mapa de calor por distrito e ano
ggplot(media_anual_distrito, aes(x = Ano, y = reorder(Distrito, -Media_Taxa_Hipertensao_Controlada))) +
  geom_tile(aes(fill = Media_Taxa_Hipertensao_Controlada), color = "white") +
  scale_fill_gradient(low = "white", high = "darkgreen") +
  labs(title = "Heat Map: Controlled Hypertension Rate by District and Year",
       x = "Year",
       y = "District",
       fill = "Rate") +
  theme_minimal()

#Vamos ajustar uma regressão linear da taxa média anual vs ano para cada distrito e ver o coeficiente da inclinação: 
# Converter o ano para numérico
media_anual_distrito$Ano_Num <- as.numeric(media_anual_distrito$Ano)

# Regressão linear por distrito
tendencias <- media_anual_distrito %>%
  group_by(Distrito) %>%
  summarise(
    inclinacao = coef(lm(Media_Taxa_Hipertensao_Controlada ~ Ano_Num))[2],
    p_valor = summary(lm(Media_Taxa_Hipertensao_Controlada ~ Ano_Num))$coefficients[2, 4]
  ) %>%
  arrange(inclinacao)

print(tendencias)

# Aplica o teste de Mann-Kendall globalmente
mk.test(media_anual_distrito$Media_Taxa_Hipertensao_Controlada)

    Mann-Kendall trend test

data:  media_anual_distrito$Media_Taxa_Hipertensao_Controlada
z = 0.13429, n = 54, p-value = 0.8932
alternative hypothesis: true S is not equal to 0
sample estimates:
           S         varS          tau 
1.900000e+01 1.796700e+04 1.327743e-02 

modelo_global <- lm(Media_Taxa_Hipertensao_Controlada ~ Ano_Num, data = media_anual_distrito)
summary(modelo_global)

Call:
lm(formula = Media_Taxa_Hipertensao_Controlada ~ Ano_Num, data = media_anual_distrito)

Residuals:
    Min      1Q  Median      3Q     Max 
-2966.1  -935.5   -27.2  1123.1  4027.1 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -27207.58  431123.04  -0.063    0.950
Ano_Num         17.59     213.11   0.083    0.935

Residual standard error: 1279 on 52 degrees of freedom
Multiple R-squared:  0.000131,  Adjusted R-squared:  -0.0191 
F-statistic: 0.006813 on 1 and 52 DF,  p-value: 0.9345

## COMPARAÇÃO COM AS TAXAS NOS ANOS COMUNS

# AVC
avc_regiao <- avc %>%
  group_by(`distrito_da_ocorrencia`) %>%
  summarise(Total_AVC = sum(as.numeric(`n_de_registos_via_verde_avc`), na.rm = TRUE)) %>%
  rename(Distrito = `distrito_da_ocorrencia`)

# Hipertensão
hipertensao_total <- hipertensao_tempo %>%
  group_by(Distrito) %>%
  summarise(Total_Hipertensao_Controlada = sum(Total_Hipertensao_Controlada, na.rm = TRUE)) 

# Juntar
df_relacao <- full_join(avc_regiao, hipertensao_total, by = "Distrito") %>%
  mutate(
    Populacao = populacao_distritos[Distrito],
    Taxa_AVC = (Total_AVC / Populacao) * 100000,
    Taxa_Hipertensao_Controlada = (Total_Hipertensao_Controlada / Populacao) * 100000
  ) %>%
  drop_na()
ggplot(df_relacao, aes(x = Taxa_Hipertensao_Controlada, y = Taxa_AVC)) +
  geom_point(color = "darkblue", size = 3) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  geom_text(aes(label = Distrito), vjust = -0.8, size = 3) +
  labs(title = "Relationship between Controlled Hypertension Rate and Stroke Rate by District",
       x = "Controlled Hypertension Rate per 100,000 inhabitants",
       y = "Stroke Rate per 100,000 inhabitants") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

# ver se tenho outliers 

boxplot(df_relacao$Taxa_Hipertensao_Controlada,
        main = "Boxplot of Controlled Hypertension Rate by District",
        ylab = "Rate")

# Boxplot for Stroke Rate
boxplot(df_relacao$Taxa_AVC,
        main = "Boxplot of Stroke Rate by District",
        ylab = "Rate")

# Para Taxa_Hipertensao_Controlada
x <- df_relacao$Taxa_Hipertensao_Controlada
Q1 <- quantile(x, 0.25, na.rm = TRUE)
Q3 <- quantile(x, 0.75, na.rm = TRUE)
IQR <- Q3 - Q1

# Limites inferior e superior
lim_inf <- Q1 - 1.5 * IQR
lim_sup <- Q3 + 1.5 * IQR

# Ver quais são os outliers
outliers <- x[x < lim_inf | x > lim_sup]
print(outliers)

named numeric(0)

#Não tem outliers
#
# Correlação
cor.test(df_relacao$Taxa_Hipertensao_Controlada, df_relacao$Taxa_AVC, method = "pearson")

    Pearson's product-moment correlation

data:  df_relacao$Taxa_Hipertensao_Controlada and df_relacao$Taxa_AVC
t = -1.4182, df = 16, p-value = 0.1753
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6929269  0.1572414
sample estimates:
       cor 
-0.3341598 

# Regressão
summary(lm(Taxa_AVC ~ Taxa_Hipertensao_Controlada, data = df_relacao))

Call:
lm(formula = Taxa_AVC ~ Taxa_Hipertensao_Controlada, data = df_relacao)

Residuals:
   Min     1Q Median     3Q    Max 
-68.30 -28.38 -18.96  27.39 106.73 

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)   
(Intercept)                  3.016e+02  9.846e+01   3.063  0.00743 **
Taxa_Hipertensao_Controlada -6.371e-04  4.493e-04  -1.418  0.17533   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 49.4 on 16 degrees of freedom
Multiple R-squared:  0.1117,    Adjusted R-squared:  0.05614 
F-statistic: 2.011 on 1 and 16 DF,  p-value: 0.1753

# Juntar por Distrito e Mês


avc_regiao_tempo <- avc_regiao_tempo %>%
  rename(Distrito = distrito_da_ocorrencia)


dados_mensais_combinados <- inner_join(avc_regiao_tempo, hipertensao_tempo,
                                       by = c("Mes_Ano", "Distrito"),
                                       suffix = c("_AVC", "_Hipertensao")) %>%
  select(Mes_Ano, Distrito, Taxa_AVC, Taxa_Hipertensao_Controlada)

# Correlação de Pearson por distrito
correlacoes_distrito <- dados_mensais_combinados %>%
  group_by(Distrito) %>%
  summarise(
    correlacao = cor(Taxa_AVC, Taxa_Hipertensao_Controlada, use = "complete.obs", method = "pearson"),
    .groups = "drop"
  ) %>%
  arrange(correlacao)

print(correlacoes_distrito)

ggplot(correlacoes_distrito, aes(x = reorder(Distrito, correlacao), y = correlacao)) +
  geom_col(fill = "purple") +
  geom_text(aes(label = round(correlacao, 2)), hjust = -0.2, size = 3) +
  coord_flip() +
  labs(title = "Correlação entre Taxa de AVC e Hipertensão Controlada por Distrito",
       x = "Distrito", y = "Correlação (r)") +
  theme_minimal()

# Juntar os dados mensais
dados_mensais_combinados <- inner_join(avc_regiao_tempo, hipertensao_tempo,
                                       by = c("Mes_Ano", "Distrito"),
                                       suffix = c("_AVC", "_Hipertensao")) %>%
  select(Mes_Ano, Distrito, Taxa_AVC, Taxa_Hipertensao_Controlada)

# Calcular correlação e p-valor por distrito
correlacoes_distrito <- dados_mensais_combinados %>%
  group_by(Distrito) %>%
  summarise(
    r = cor(Taxa_AVC, Taxa_Hipertensao_Controlada, use = "complete.obs", method = "pearson"),
    p_value = cor.test(Taxa_AVC, Taxa_Hipertensao_Controlada)$p.value,
    .groups = "drop"
  ) %>%
  arrange(r)

# Ver distritos com correlação negativa significativa
correlacoes_significativas <- correlacoes_distrito %>%
  filter(r < -0.3, p_value < 0.05)

print(correlacoes_significativas)

ggplot(correlacoes_distrito, aes(x = reorder(Distrito, r), y = r, fill = p_value < 0.05 & r < -0.3)) +
  geom_col() +
  geom_text(aes(label = round(r, 2)), hjust = -0.2, size = 3) +
  coord_flip() +
  scale_fill_manual(values = c("grey70", "darkred"), labels = c("Não significativa", "Significativa"), name = "Correlação negativa significativa") +
  labs(title = "Correlação entre Hipertensão Controlada e AVC por Distrito",
       x = "Distrito", y = "Coeficiente de Correlação (r)") +
  theme_minimal()

NA
NA

# Correlação entre hipertenção e avc com as taxas
# Considerando o maior período de tempo dado pela data

#população média pr distritos entre 2015 e 2023
populacao_media <- c(
  Lisboa = 2276293,
  Porto = 1795549,
  Setubal = 869770,
  Braga = 843697,
  Aveiro = 703022,
  Faro = 462522,
  Leiria = 457242,
  Santarem = 437520,
  Coimbra = 411935,
  Viseu = 362275,
  Madeira = 251025,
  Acores = 238127,
  Viana_do_Castelo = 234672,
  Vila_Real = 195422,
  Castelo_Branco = 183996,
  Evora = 155441,
  Beja = 144772,
  Guarda = 148932,
  Braganca = 127356,
  Portalegre = 109975
)

hipertensao_tempo_all <- hyp_final %>%
  mutate(Mes_Ano = format(ym(tempo), "%Y-%m")) %>%
  group_by(Mes_Ano, distritos.x) %>%
  summarise(
    Total_Hipertensao_Controlada = sum(as.numeric(total), na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    Populacao = populacao_media[distritos.x],
    Taxa_Hipertensao_Controlada = (Total_Hipertensao_Controlada / populacao_media) * 100000
  )

Warning: There was 1 warning in `stopifnot()`.
ℹ In argument: `Taxa_Hipertensao_Controlada = (Total_Hipertensao_Controlada/populacao_media) *
  1e+05`.
Caused by warning in `Total_Hipertensao_Controlada / populacao_media`:
! longer object length is not a multiple of shorter object length

hipertensao_total_2015 <- hipertensao_tempo_all %>%
  group_by(distritos.x) %>%
  summarise(Taxa_Hipertensao_Controlada = sum(Taxa_Hipertensao_Controlada, na.rm = TRUE)) %>%
  arrange(desc(Taxa_Hipertensao_Controlada))%>%
  rename(Distrito = distritos.x)


avc_hipertensao_all <- avc_regiao_tempo %>%
  inner_join(hipertensao_total_2015, by = "Distrito")

ggplot(avc_hipertensao_all, aes(x = Taxa_Hipertensao_Controlada, y = Taxa_AVC)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Correlation between Controlled Hypertension Rate and Stroke Rate by District",
       x = "Controlled Hypertension Rate from 2015 to 2024",
       y = "Stroke Rate from 2022 to 2024") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

correlacao <- cor(avc_hipertensao_all$Taxa_Hipertensao_Controlada, avc_hipertensao_all$Taxa_AVC)
print(correlacao)

[1] -0.1876239

#Testes não paramétricos
#correlação de spearman - monótona não linear
c <- cor(avc_hipertensao_all$Taxa_Hipertensao_Controlada, avc_hipertensao_all$Taxa_AVC, method = "spearman")
print(c)

[1] -0.1545749

# medida não paramétrica que avalia a associação monotônica entre duas variáveis. Ou seja, verifica se, à medida que uma variável aumenta, a outra tende a aumentar (ou diminuir), sem exigir uma relação linear.
#existe uma correlação neagtiva fraca (-0.2125704) entre as duas variáveis


#coreelação de kendall - mede a força e a direção da relação monotônica
k <- cor(avc_hipertensao_all$Taxa_Hipertensao_Controlada, avc_hipertensao_all$Taxa_AVC, method = "kendall")
print(k)

[1] -0.1179419

resultado <- cor.test(avc_hipertensao_all$Taxa_Hipertensao_Controlada, avc_hipertensao_all$Taxa_AVC,
                      method = "kendall")
print(resultado)

    Kendall's rank correlation tau

data:  avc_hipertensao_all$Taxa_Hipertensao_Controlada and avc_hipertensao_all$Taxa_AVC
z = -3.7095, p-value = 0.0002077
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
-0.1179419 

#existe uma correlação negativa fraca (-0.1589298) entre as duas variáveis
# 1.72e-06 < 0,05: A correlação é estatisticamente significativa, embora fraca.

#correlação AVC e hipertensão SÓ nos anos em comum


avc_hipertensao <- avc_regiao  %>%
  inner_join(hipertensao_total, by = "Distrito")

ggplot(avc_hipertensao, aes(x = Total_Hipertensao_Controlada, y = Total_AVC)) +
  geom_point(color = "steelblue", size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Correlação entre Hipertensão Controlada e Registos de AVC por Distrito",
       x = "Total de Hipertensos Controlados",
       y = "Total de Registos de AVC") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

correlacao <- cor(avc_hipertensao$Total_Hipertensao_Controlada, avc_hipertensao$Total_AVC)
print(correlacao)

[1] 0.9733902

```